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Abstract. We investigate the dynamics of capillary filling using two lattice Boltzmann 
schemes: a liquid-gas model and a binary model. The simulation results are compared to the 
well-known Washburn's law, which predicts that the filled length of the capillary scales with 
time as I oc t 1 ^ 2 . We find that the liquid-gas model does not reproduce Washburn's law due 
to condensation of the gas phase at the interface, which causes the asymptotic behaviour 
of the capillary penetration to be faster than t 1 ^ 2 . The binary model, on the other hand, 
captures the correct scaling behaviour when the viscosity ratio between the two phases is 
sufficiently high. 



1 Introduction 



The aim of this paper is to present an effective modelling method for simulating liquid which fills a 
capillary tube that initially contains gas. Capillary tubes have a hydrophilic inner surface and the energy 
liberated in wetting this surface is used to drive the fluid up the tube. 
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Fig. 1. The two simulation setups, (a) A drop resting at equilibrium on a surface. (8 is the contact angle.) (b) 
Capillary flow along a tube. I is the length of tube filled with fluid and l eii — I + H/2 is an effective filled length. 
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The classical analysis of capillary penetration is due to Washburn (T). He assumed that the liquid is 
incompressible and that the fluid flow has a parabolic profile (Poiseuille flow). In two dimensions, the 
average velocity of a parabolic flow is 

H 2 dp 

v = — (1) 

12?7 dx 

where H is the capillary tube width, 77 is the liquid viscosity, and dp/dx is the pressure gradient that sets 
up the flow. The Laplace pressure drop across a curved interface of radius R is given by Ap = jig/R, 
where ji g is the liquid-gas surface tension. The pressure gradient in the fluid is therefore 

dp = _1lg (2) 
dx Rl 

where I is the length of the liquid column that has penetrated the capillary. R is related to the dynamic 
contact angle 9 (see Fig.QJb) for a definition of 6) through R = H/2 cos 9. By substituting Eq. © into 
(fl~|l and using v = dl/dt we obtain 

,_(3te|p») 4 (t + 4o ,l (3) 

where to is an integration constant. An alternative way to derive Eq. (O is to equate the dissipation of 
energy due to viscosity to the energy liberated as the liquid wets the surface. In the original analysis, 
Washburn neglected the viscous dissipation of energy in the gas phase (the viscosity ratio between water 
and air is around ~ 10 3 so this is a very good approximation in this case) and the deviation from the 
Poiseuille flow velocity profile at the inlet and near the curved interface 12. 

In this paper we consider two approaches to modelling capillary dynamics. The first is a van der 
Waals liquid-gas model and the second is a binary fluid model with a viscosity difference between the 
two phases. We show that the liquid-gas model does not reproduce Washburn's law due to condensation 
of the gas phase at the interface, which causes the asymptotic behaviour of the capillary penetration to be 
faster than t 1 / 2 . The binary model, however, captures the correct scaling behaviour when the viscosity 
ratio between the two phases is sufficiently high. Other authors [3 4 5 6] have discussed different lattice 
Boltzmann approaches to model capillary filling. Their results are broadly similar to those reported here, 
but the models differ in the details of the behaviour of the dynamic contact line. 



2 A liquid-gas model 

The pressure tensor for a liquid-gas system, resulting from a Landau free energy functional, is 1718 



where 



P a p = [po - np\7 2 p- ^\Vp\ 2 ) S a /3 + Kd a pdpp : (4) 



Po = ap 2 (5) 

1 — bp 

is the van der Waals bulk pressure, p is the fluid density and k is a parameter related to the surface 
tension. This leads to liquid-gas phase separation below a critical temperature. For this investigation 
we chose the parameters k = 0.02, a = 9/49, b = 2/21 and T = 0.56 giving a surface tension of 
jig = 0.0112 and liquid and gas densities of pi = 4.54 and p g = 2.59, respectively. 



2.1 Lattice Boltzmann implementation 

We used a two-dimensional, nine velocity vector, free energy, multiple-relaxation-timescale lattice 
Boltzmann method [9]. Here we briefly outline the method and refer the reader to ifTol for more de- 
tails. The system is divided up into a square grid of nodes, and on each node there is a particle dis- 
tribution function /j(r, t). The label i denotes a particular lattice velocity ej, defined by eo = (0, 0), 
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Gradient in p at the Boundary Gradient in § at the Boundary 



Fig. 2. The equilibrium contact angle as a function of the equilibrium gradient in p or <j> at the boundary for (a) the 
liquid-gas system and (b) the binary system. The circles were obtained numerically by quasistatically scanning the 
equilibrium gradient in time, using v = 1/6 in (a) and vi — 0.83 and v g = 0.067 in (b). The dashed curve comes 
from Young's law <[8l», based on numerically calculated surface tensions. 

ei.2 = (±c, 0), 63,4 = (0, ±c), eg^ = (±c, ±c), and e7 g = (T c > ± c )- The lattice speed c is given by 

„ _ Ax 
c ~ Af 

The time evolution equation for the particle distribution function is 

f (r + eAt, t + At) = f (r, t) - M _1 SM [f - i eq ) , (6) 

where /,; has been written as a column vector, M is a matrix that performs a change of basis, S is a 
diagonal matrix which defines different relaxation times for different modes and f eq is an equilibrium 
distribution function. Details of a suitable choice for M, S and f eq are given in AppendixlAl 

In the limit of long length and timescales, Eq. © leads to the continuum Navier-Stokes equation 

d t {pv a ) + dp(pv a vp) = -dpP a p + dp (up [dpv a + d a vp]) , (7) 

which determines the dynamics of the system. The parameter u, the kinematic viscosity, is related to the 
dynamic viscosity by -q = pu. 



2.2 Numerical results 

To check the equilibrium properties of the model, we numerically measured the contact angle and com- 
pared it with theory. Figure[Tfa) shows a drop resting on top of a solid surface. In general, the solid-liquid 
jsi, solid-gas 7 SS and liquid-gas -fi g surface tensions are different. At the contact point A the balance of 
forces gives Young's law: 

cosQeq = leg - Isl 
llg 

where 9 eq defines the equilibrium contact angle. 

We used a system of size 300 x 50 lattice units and placed an initially semi-circular drop on the 
lower wall. Non-slip boundary conditions were simulated by using a bounce-back scheme as well as 
setting boundary nodes to rest after each streaming step (this and the MRT-LB scheme were found to be 
necessary to stop unphysical currents appearing near to the interfaces [ 10|). The wetting properties of the 
surface were changed by altering the gradient in p, d n p\ b , as it appears in the equilibrium distribution, 
at the boundary [11]. Figure|2a) shows that there is good agreement between the numerically measured 
contact angle (the circles) and that predicted by Young's law (the dashed line). 

Next, we focus on simulating a dynamical system, that of capillary filling, and test the Washburn 
relation in Eq. (01. The simulation setup is illustrated in Fig.QJb). The system consists of a lattice of size 
700 x 40 lattice units with periodic boundary conditions in the x direction. The upper and lower sides of 
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Fig. 3. The distance of the liquid-gas interface along the capillary as a function of time. Circles are simulation 
results using v = 1/6, the solid line is Washburn's law (O and the dashed line is Washburn's law taking into 
account dissipation in the gas phase. 
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Fig. 4. (a) The velocity field around the liquid-gas interface during capillary filling using the standard liquid-gas 
model, (b) The x component of the velocity as a function of x at the centre of the capillary. 



the system have two sets of boundary conditions. In the middle portion the boundaries are non-slip and 
wetting (denoted by the hashed region in the diagram) and this represents the sides of the capillary. The 
length of this is set to be L — 350. At either side of the capillary the boundary conditions are periodic 
in the y direction, and hence allow slip, and these areas represent a reservoir of liquid and gas. 

Using the results from Fig. [2ja) we chose the equilibrium gradient at the boundary to be d n p\ b = 
—0.144 such that the equilibrium contact angle is 9 eq — 60°. This wetting interaction induces the liquid 
and gas interface to form a meniscus. The system is initialised such that I = 10 and evolved in time 
using the lattice Boltzmann algorithm. The solid line in Fig. [3] shows a plot of I against time. When 
compared against the theoretical prediction of Washburn (given by the dashed line) we observe that 
agreement is poor. One reason behind this might be that viscous dissipation in the gas phase cannot be 
ignored. However, taking this into account (as shown be the dotted line) makes the predicted flow even 
slower. 

The reason behind the discrepancy between theory and simulation can be clearly seen if we observe 
the fluid flow profile near to the meniscus, as plotted in Fig. Ufa). We find that whilst fluid in the liquid 
phase is moving to the right and filling up the capillary, the fluid in the gas phase is moving left and 
condensing to form liquid, which helps to significantly increase the speed of the interface. This is seen 
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even more clearly in Fig. 12b), which shows the x component of the fluid velocity as a function of 
distance down the tube. It is because the system primarily exhibits a condensation driven interface 
velocity, rather than sucking fluid through the capillary, that Washburn's law breaks down. 

3 A binary model 

3.1 Lattice Boltzmann implementation 

One way to prevent condensation is to simulate the system as a binary fluid B7I8L We associate A 
particles with the liquid phase and B particles with the gas phase. Because particle species is strictly 
conserved, condensation is no longer permitted. Obtaining a viscosity ratio between the two phases is 
achieved by making the kinematic viscosity v a function of the order parameter <fr = (pa — pb)/{pa + 
PbY 

v = v g + i ^r(yi- v g ) (9) 

such that the viscosity has the bulk values v\ and v g in the liquid and gas phases, respectively. 

We choose the well known "0 4 theory" to model the phase separation into two distinct phases 
(f> = ±1. This has a pressure tensor and chemical potential given by 

P af3 = (pa - n4>V 2 4> - ||V</)| 2 ) S aP + Kd a <j>d p <t>, (10) 

ft= A(-0 + 3 ) -kV 2 0. (11) 
The bulk pressure in this case is 

p Q = p^ + A(-tj 2 + ^ 4 ). (12) 

In this study, we use the parameters A = 0.04 and k — 0.04, which lead to an interface width of ~ 4 
lattice sites and a liquid-gas surface tension of ji g = 0.0389. The average density of the system was 
taken to be p = 1. 

As well as a time evolution equation for /j in Eq. which gives the Navier-Stokes equation, there 
is a second lattice Boltzmann equation 

S {r + eAt,t + At) = g e9 (r,i) (13) 

which describes the time evolution of gi(r, t), the order parameter distribution function. An expression 
for the equilibrium distribution g eq is given in Appendix|A] This leads to an advection diffusion equation 
for the order parameter 

^+9 Q (K)=MV 2 /i, (14) 

where M is a mobility parameter. 

3.2 Numerical results 

Firstly, the equilibrium properties of the system were verified lfl2l . Figure |2b) shows the contact angle 
obtained numerically as a function of the equilibrium gradient of <fi at the boundary. As with the liquid- 
gas model, good agreement is found with Young's law. (Note that the correct curve, in the case when 
there is a viscosity difference between the liquid and gas phases, is not produced if a standard BGK 
lattice Boltzmann algorithm is used. Details of why this is the case are presented in [ 10|.) 

Figure|5ja) shows a plot of capillary filling distance I as a function of time using the binary model. 
The circles are results when the liquid and gas phases have equal viscosity. They form a straight line 
because viscous dissipation occurs at approximately the same rate at any given point down the tube, 
and so the total dissipation is independent of the position of the interface. The solid, straight line shows 
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Fig. 5. (a) The distance of the liquid-gas (A-B) interface along the capillary as a function of time for the binary 
fluid model. Circles: vi — v g = 1/6, squares: vi = 1.17 and v g = 0.017, and the solid and dashed lines are 
theoretical predictions (see text for details), (b) The corresponding variation in dynamic contact angle. The dashed 
line extrapolates the results to the limit —> oo. 
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Fig. 6. (a) The velocity field around the liquid-gas (A-B) interface during capillary filling using the binary model 
{vi = 0.83, v g = 0.067). (b) The x component of velocity as a function of x along the centre of the capillary. Solid 
line: v% = v g = 1/6, dashed line: vi = 0.83, v g — 0.067. 



the theoretically expected profile based on Poiseuille flow down a tube of length L = 350 using the 
numerically measured contact angle of 9 = 63° (see below). The agreement is not exact as we have 
ignored dissipation at the inlet and outlet of the tube. This can be taken into account by extending the 
effective length of the tube either end by ~ H/2. Subsequently, we use an effective length of filling 
l & >' = I + H/2, as shown in Fig.[TJb). The dashed, theoretical curve in Fig.Qa) takes into account this 
refinement and shows much closer agreement with the simulation results. 

Measurements of the meniscus contact angle 9 (see Fig. |TJb)) were made by performing a least 
squares fit of the interface profile to a circular section. The circles in Fig. |5jb) show 6 as a function 
of Hjl e *'. (Note that increases in time so this plot starts on the right side of the diagram and 
moves left.) The noise in this curve is not real but a result of the measurement technique. We observed 
that 9 rapidly reaches a stable value of 9 ~ 63°. This is significantly different from the equilibrium 
contact angle of 9 eq — 60°. The reason behind this is revealed if we look at the flow field near to the 
meniscus, as plotted by the black arrows in Fig. |6ja). In the bulk phases, the flow down the capillary 
tube is parabolic. On the other hand, the interface itself moves at a constant velocity, except very close 
to the boundaries where the non-slip boundary conditions prevent this (in this case diffusion allows the 
"slip"). The transition between these two profiles is achieved by fluid in the liquid phase being driven 
into the corners, towards the contact points, and, conversely, the gas phase being pushed away (observe 
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Fig. 7. (a) The variation in the velocity of capillary filling as a function of mobility M, using v\ = v g = 1 /6. (b) 
The corresponding change in dynamic contact angle. 

the arrows near to the contact points in Fig. HI a)). The driving force for this process is the difference 
between the dynamic and equilibrium contact angles A9 — 9 — 9 eq . Another way of looking at this is 
to say that, because energy is dissipated in the flow field near the interface, the dynamic contact angle 
9 that appears in the Washburn relation (01 must be greater than 9 eq to reflect the fact that there is less 
available energy to do useful work (in this context useful work means the energy for pulling the fluid 
into the capillary). 

The curves in Fig. |6jb) clearly show that the velocity in the liquid and gas phases are identical, as 
compared to the very different situation that was observed in the liquid-gas system in Fig. Hfb). The 
dip in velocity at the interface is because of the change in flow profile, as discussed in the previous 
paragraph. The dashed curve is for a high viscosity ratio, and the bump in the velocity profile at the end 
of the capillary resulted from the formation of vortices indicative of a moderate Reynolds number in the 
gas phase. 

The squares in Fig. |5ja) show the results of a simulation with a large viscosity ratio between the 
two phases of vi/v g = 70, sufficiently high that viscous dissipation in the gas phase can, justifiably, be 
ignored. The dashed curve is from Washburn's law based on a dynamic contact angle of 9 = 63°. Very 
good agreement is achieved at late times. The difference at early times can be explained if we examine 
the dynamic contact angle in Fig. [2b). In this case, A9 starts off large and gradually decreases as l e '$ 
gets bigger. The value 9 — 63° was used above because it is representative of the late stages of the 
simulation. If we set 9 — 68° in Washburn's law then we obtain the dotted curve, which can be fitted 
well to the early time behaviour. The dashed line in Fig. |5Jb) is a linear extrapolation of the contact 
angle data showing that, as the liquid penetrates far into the tube and the meniscus velocity tends to 
zero, 9 = 9 eq as expected. 



3.3 The role of particle mobility, M 

Unlike the liquid-gas model, the binary model has an extra parameter, M, which determines the diffusive 
behaviour of the liquid (A) and gas (B) particles. To assess the effect of varying M on the system, 
we consider capillary filling for the case when the the viscosity of the two phases is the same. This 
particularly simple case was chosen because the velocity of filling is constant and can be calculated 
from a linear fit to the filling profile. Figure [3 a) shows the velocity of filling against M. In the limits 
of small or high M the algorithm was found to be numerically unstable. In the stable region there is a 
general trend of increasing velocity with increasing M. This is rather an unexpected result. The right 
hand side of the diffusion equation (TBI) is normally associated with dissipation of energy towards a 
free energy minimum ([Pil l and yet we find that increasing its size increases the energy available to do 
work (in this case work in filling the capillary tube). The reason that the contact angle decreases with 
increasing M is that the diffusion alleviates the no-slip boundary conditions - for larger M the interface 
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slips more easily across the surface. Thus dissipation due to flow in the vicinity of the meniscus, and 
hence AO, are reduced. 



4 Summary 

We have assessed two different ways of simulating capillary flow. Firstly, using a liquid-gas approach, 
the filled length of a capillary tube was found to increase much more rapidly than predicted by Wash- 
burn's law. This discrepancy was due to condensation of the gas phase at the interface. This is allowed 
within the formalism we used because the van der Waals equation of state describes a liquid in equi- 
librium with its vapour. Moreover, the effect of condensation is large because the liquid and gas have 
similar densities. Essentially we are modelling a system close to its critical point. Results in [5| show 
that a good fit to the Washburn equation is obtained if the liquid-gas density ratio is large. 

We then considered a binary fluid, comprising two different types of particles, where evaporation 
and condensation is not permitted. Very good agreement was achieved with the theoretical expression 
in Eq. (O, as long as the dynamic contact angle was used in the fit. We argued that the dynamic contact 
angle differed from the equilibrium value because of dissipation of energy from the flow field near to 
the meniscus. 

We hope that this simulation technique will prove useful in studying capillary filling in porous media 
and in microchannels with differing geometries or surface structures and patterning. 
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A Details of the lattice Boltzmann scheme 

This appendix gives details of the lattice Boltzmann scheme used for this study. The matrix M, that 
describes a change of basis, is given by 
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(15) 



In the new basis the diagonal matrix 

S = diag (0,1, 1,0, 1,0,1, (16) 

sets the relaxation rates of different modes. Some of these are arbitrarily set to zero and these correspond 
to conserved quantities, e.g. the top line dotted with f gives the density p. The quantity u> = 2/ (Qv + 1) 
sets the kinematic viscosity. Because v is a function of <f> in Eq. (0 it might seem necessary to calculate 
the collision matrix, M _1 SM, at each node at each time-step, which would be computationally very 
slow. Our approach is to make a look up table containing ~ 10 4 matrices with different values of 
viscosity and simply pick the closest match. 
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The equilibrium distribution can be written in the form 

ft q i v ) =^t\Po- npV 2 p + e ia pu a + 

(pu a Uf3 + A [u a dpp + updap + 5 al3 u 1 d 1 p}) 
+ §t (wf x d x pd x p + wfdypdyp + w xy d x pd y p 



Saf 



(17) 



for i = 1, .., 8, where W1.4 — |, u>5_§ = j^, and summation over repeated indices is assumed. Other 



parameters are wf 2 
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and Wyg = 
The i 



— i. This choice was made to reduce spurious velocities generated at interfaces lfl3ll . 
stationary value is chosen to conserve mass: 



»=i 



For the binary model, the equilibrium distribution is given by 

8 



(18) 



(jm a U/3 



(19) 



For the order parameter distribution functions, the relaxation rates are all set to 1 . During the lattice 
Boltzmann procedure, it is necessary to numerically calculate both derivatives (e.g. d x p in the equilib- 
rium distribution ( fTTI i) and the Laplacian (e.g. to obtain the chemical potential (fTTT)). These continuous 
quantities are calculated from stencils, discrete operators which use neighbouring lattice sites. The best 
choice of stencils to reduce spurious velocities is given by: 
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(20) 
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